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We present a generic reweighting method for nonequilibrium Markov processes. With nonequilib- 
rium Monte Carlo simulations at a single temperature, one calculates the time evolution of physical 
quantities at different temperatures, which greatly saves the computational time. Using the dy- 
namical finite-size scaling analysis for the nonequilibrium relaxation, one can study the dynamical 
properties of phase transitions together with the equilibrium ones. We demonstrate the procedure 
for the Ising model with the Metropolis algorithm, but the present formalism is general and can be 
applied to a variety of systems as well as with different Monte Carlo update schemes. 
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The Monte Carlo simulation has served as a stan- 
dard method to treat many body problems in statistical 
physics Metropolis et al. used the importance 
sampling, which generates the states with a probabil- 
ity proportional to the Boltzmann factor. The impor- 
tance sampling method can be considered as a reweight- 
ing technique in a sense that one samples a trial distribu- 
tion and takes an average for the true distribution with 
reweighting. A more direct way of reweighting is found 
in the histogram reweighting method due to Ferrenberg 
and Swendsen fsj , where equilibrium thermal averages for 
a range of temperatures can be calculated from a single 
simulation. This method greatly improved the efhcien- 
cies of Monte Carlo simulations, however, the histogram 
reweighting technique has to be used only to calculate 
equilibrium properties. 

Recently, nonequilibrium relaxation method has been 
successfully applied to the study of critical phenom- 
ena U, |5|, ly, . In the nonequilibrium relaxation method, 
simulations were performed for several temperatures; the 
critical temperature, the dynamical exponent and other 
quantities are estimated using the scaling behavior of 
nonequilibrium process. If we combine the strength of 
nonequilibrium relaxation method with a reweighting 
technique, we can expect an effective method of simu- 
lation. 

In this paper, we will present a generic reweight- 
ing method that is applicable to both equilibrium and 
nonequilibrium systems. With reweighting at nonequi- 
librium, only a simulation at a single temperature is re- 
quired. Our method of reweighting may also be appli- 
cable in a multitude of other nonequilibrium processes, 
such as the contact process [sj and the driven diffusive 
systems 

Reweighting at nonequilibrium can be easily under- 
stood as follows. Consider every simulation up to the tih 
Monte Carlo step as a sequence of states (or path). 



(0-1,0-2, • • • , CTt) 



(1) 



where at is the system configuration at time t. Hereafter, 
we refer to Monte Carlo step simply as the time of sim- 
ulation. Such path Xt can be generated using any Monte 



Carlo method at a temperature T . The objective is to 
calculate the relative probability of generating the same 
path Xt at another temperature T' . 

Suppose many simulations were performed at an in- 
verse temperature /3 — l/k^T to obtain a set of paths 
xItJ = li'-'i'^- (From now on, (3 will be sometimes 
referred to as temperature). Dynamical thermal aver- 
age of some quantity Q{t) can be obtained by {Q{t))i3 — 
X]j=i Q(^t )• To calculate the dynamical thermal 
average at another temperature /?', the measured quan- 
tity Q{xl) has to be reweighted with a set of weights wl. 
For the same set of paths , the thermal average at j3' 
is obtained as 



(2) 



Although not labeled explicitly in Eq. the set of 
weights wl depend on the simulation temperature (3 
and the new temperature /?'. To calculate the weights, 
the following algorithm is implemented for each path 
^1,3 = l,•■■,'^■ 
l. Assume that the Monte Carlo simulation is car- 
ried out at a temperature (3, and a path xl = 
{a^, • ■ • , Oj ) for some arbitrary time t is obtained. 

2. To go from time t, let a'^ be a trial configuration 
and T(o-'"' \al) be the probability to select this con- 
figuration. If the trial move is accepted with the 
acceptance probability A^(o-'"'|o-(), then the new 
configuration aX t + \ becomes CTj^j^ — a'^ . If the 



move is not accepted, o-(_|_j^ 



3. In terms of the transition probability 



T{a''' \a-j.)Ap{a'-' \al) if the move is accepted and, 
Pp{<jl+,\<ji) = T{a''\al)[l-A0{a''\ai)] otherwise. 

4. Then the appropriate weights wl^i can be obtained 



we may write as Pf3{cr-l_^_^\af) 



2 



by 



w: 



t+i 



7^ —'t 



with 



(3) 



5. Repeat steps 2 to 4 until t reaches some predeter- 
mined maximum simulation time. 

For each path xl,j = 1, ■ • ■ ,n, these steps are repeated. 
Step 2 is simply the standard Monte Carlo procedure of 
acceptance and rejection. Step 3 calculates the transition 
probability after the decision on acceptance was made. 
Step 4 updates the appropriate weight. In this method, 
essentially the usual Monte Carlo update is carried out, 
and at the same time the weights were updated. 

The proof of Eq. involves two ingredients, the Se- 
quential Importance Sampling (SIS) [lOllllLll^ method 
and a generic Monte Carlo method. A brief description of 
SIS method will be discussed here. The SIS method as- 
sumes that the distribution for a multi-dimensional state 
Xt — (<7i , o'2 , • ' ' J '''t ) is known and let it be denoted by 
T:t{xt)- Implementation of the SIS algorithm is as follows: 

1. Draw cF^^i from a trial distribution gt+i{cr-l^-^\xt) 



to form 



2. The importance weight at t -I- 1 is given by 



7rt+i(x^+i) 



7rt(S't).9t+i(o-t+il^t) 



(4) 



The key to understanding the reweighting method is to 
realize that in a Monte Carlo simulation, the time se- 
quence xl is sampled from the true probability distribu- 
tion of the path at some temperature (3. The Monte Carlo 
method is Markovian and the probability of obtaining the 
path is given by 

Pf,{xi) = Pf3iaiWl_,) ■ ■■Pp{ai\a{)P{ai) (5) 

where Ppicl {(^i-i) ^^e probability of getting a system 
configuration al at time t given that the system configu- 
ration is at an earlier time t—1. P{<t{) is the prob- 
ability of choosing the initial configuration. One may 
obtain the probability distribution function, P^' (xj ) , of 
the path at another temperature /3' using Ppix^) as a 
trial distribution. The SIS method can then be used by 
defining the following quantities in Eq. Q as 



TTtixl) = Pp.{xl) 

gt+i{(Ti+i\xl) EE Pp{alj^^\xl) 



(6) 
(7) 
(8) 



The requirement for the Monte Carlo process being 
Markovian implies P0{alj^i\xl) = P^(cr^_|_j^ Icr^ ), and us- 
ing Eq. ||SJ), we obtain 



Pp,{xl)Pp{ai^,\xl) 



1,3 — 



Ppi^l+iWi) 



(9) 



which is as presented in Eq. (PJ. The initial condition in 



Eq. ®is' 



1 because the initial system configuration 



is chosen from a distribution independent of temperature. 

To illustrate how the reweighting algorithm is imple- 
mented, we use the ferromagnetic Ising model on a square 
lattice. Its Hamiltonian is given by 



(ki) 



(10) 



where the sum is over nearest neighbors and Sk takes the 
values ±1. Periodic boundary conditions are used on a 
L X L lattice. We use the single spin fiip update with the 
Metropolis acceptance rate, but other update schemes 
are also applicable. The simulation is performed at /3; the 
objective is to reweight it to (3' . The system configuration 
a is denoted by the states of all spins, a — {si, • • • , sjv}, 
where N is the total number of sites. The initial system 
configuration at time t = 1 is set to Sfc = -1-1 for all 
k — 1,---,N. Hence wl = 1 for j = l,---,n; recall 
that j is used for indexing different paths. For each path 
xl = ia{,ai,---:CFi), 

1. Choose a spin and flip it with probability, 
Ap{a'^\ai) = min[l,exp(-/3Ai;)]. Here, AE is the 
energy change due to the spin flip. 

2. At this point, there are three possible outcomes, 

(a) If AE' < 0, the proposed spin flip is always 
accepted; and we obtain urj.^i = wf x 1. 

(b) If AE > and the proposed spin flip is ac- 
cepted, we obtain Wj+i = wl x exp(— (/?' — 
(3)AE). 

(c) If AE > and the proposed spin fiip 
is rejected, we obtain w^^i — wl x 
(1 - exp{-f3'AE))/{l - exp{-pAE)). 

3. Although the weights should be updated for every 
single spin fiip move, thermal quantities may be 
recorded only once per TV single spin fiip steps, for 
example, to save memory. Let r be the time in unit 
of Monte Carlo steps per site (MCSS); weights w^ 

and moments of magnetization (m''^)'^,A = 1,2,4 
(i.e. m, and m'^ moments) are accumulated into 



the sums, J2j ^tj J2j (™ )r ^^'i J2j {"m'T^w: 



\\3 i 
T r • 



We carried out many simulations to obtain a set of paths 
Xj, J = 1, • • • , n. Finally, the dynamical thermal averages 
for each time r, (m(T)^)^/, for example, can be calculated 
from the accumulated sums using Eq. ((2Jl. In practice 
reweighting to many temperatures P'l, P2, ■ ■ ■ were made 
at the same time in a single run. 

Fig.fflshows plots of the temporal evolution of magne- 
tization for several reweighted temperatures from simula- 
tions at a single temperature T — 2.270 (bold line) on the 
i = 64 square lattice. These data were generated with 
three independent runs to estimate error bars; for each 
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FIG. 1: Plot of magnetization for reweighted temperatures 
from simulations at a single temperature T — 2.270 (bold 
line). From top to bottom, T =2.268, 2.269, 2.270, 2.271, 
2.272. System size is L = 64. r represents Monte Carlo time 
in units of MCSS. Insert shows the region between 1000 to 
2000 MCSS with representative error bars. 

run, averages were taken over 2.94 x 10^ samples. Effec- 
tive reweighting temperature range is about AT — 0.002 
for L = 64. Insert shows the region between 1000 and 
2000 MCSS with representative error bars. Error bars 
become larger as reweighting is done at temperatures fur- 
ther away from the simulation temperature. Similar to 
the argument put forward by Ferrenberg and Swendsen 
0, reweighting is effective when the distributions Pi3{x{) 
and Pi3'{xl) have sufficient overlaps. An additional sim- 
ulation, performed at T = 2.268, for example, yields the 
results consistent with the reweighted ones within statis- 
tical errors. A more detailed comparison of actual and 
reweighted data will be discussed shortly. 

The ratio of the moments of the order parameter is 
used for the analysis of the phase transition . Fig. |2 
shows plots of (m(T)^)/(TO(r)^)^ for several reweighted 
temperatures from simulations at T = 2.270. From 
top to bottom, T = 2.272, 2.271, 2.270, 2.269, 2,268, 
2.267. The lattice size was L = 32 and accuracy of 
reweighting was checked by performing an additional 
simulation at T = 2.268. This simulation was then 
reweighted to T = 2.270 and compared to the curve 
from the original simulation. The dashed line in Fig. 12 
shows {m{T)^) / {m{Ty)^ at T = 2.270 reweighted from 
T = 2.268. Insert shows that the difference between the 
dashed line and solid line is about 5 x lO^** while error 
bars of (TO(r)'*)/(m(T)^)^ are of the order 10"'^. Aver- 
ages were taken with 1.28 x 10^ samples. The systematic 
errors due to reweighting are smaller than the statistical 
errors for this temperature shift. 

Equilibrium properties can be calculated using 
reweighting by performing the simulation up to equilib- 
rium and beyond, and then disregard data at nonequi- 
librium. In this domain, our reweighting method yields 
the same results as the histogram reweighting Q method. 
Fig.|31shows Binder's cumulants 13] for L = 32, L = 48, 
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FIG. 2: Values of {m{T)''') / {m{T)^)^ for reweighted temper- 
atures from simulations at T = 2.270. From top to bot- 
tom, T = 2.272, 2.271, 2.270, 2.269, 2,268, 2.267. To check 
the range of reweighting, an additional simulation was per- 
formed at T = 2.268 and reweighted to T = 2.270. Dashed 
line shows {m{T)'^) / {m{Tff at T = 2.270 reweighted from 
T = 2.268. Insert shows that the difference between the 
dashed line and solid line is about 5 x 10~* while error bars 
on {m{T)*') / {mlyT)'^)'^ are of the order 10""^. The system size 
is L = 32. T represents Monte Carlo time in units of MCSS. 
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FIG. 3: Binder's cumulants for L = 32, L = 48, L = 64 and 
L = 80. Simulation temperatures are indicated by arrows, 
T = 2.27 for L = 32, 48, T = 2.269 for L = 64, 80. Error bars 
when not shown are smaller than the size of the symbols. 



i = 64 and L — SQ. Simulations were performed at one 
temperature for each lattice size. The crossing of curves 
with different sizes yields the determination of the criti- 
cal temperature Tc. In the present case, we confirm that 
the exact Tc = 2.2691 • ■ • is well reproduced. Simula- 
tion temperatures are indicated by arrows; T = 2.27 for 
L = 32, 48 and T = 2.269 for L = 64, 80. The reweighting 
temperature range we used is AT — 0.004 for L — 32, 48, 
AT = 0.002 for L = 64 and AT = 5 x 10"'^ for L = 80. 
Error bars were generated with several independent simu- 
lations, and they are smaller than the size of the symbols. 
This suggests that the effective reweighting temperature 



4 



1.20 1 

CM 




0.001 0.01 0.1 



FIG. 4: Finite-size scaling plot of (m(r)'')/(m(r)^)^ for L = 
32, 48, 64 and 80 at Tc = 2.2691 ■ ■ ■. The dynamical exponent 
is chosen as z = 2.143 for this plot. 



range is actually larger than the temperature range we 
used in our calculations. For each independent run, av- 
erages were taken with 1.28 x 10^ samples for L = 32, 
2.59 X 10^ samples for L = 48, 1.28 x 10^ samples for 
L = 64 and 4 x IC* samples for L = 80. 

The dynamical finite-size scaling can be used 
in combination with nonequilibrium relaxation [TrL IT^ . 
Making use of the scaling form pF7| . 

{m{Tf)/{m{Tff^g{TL-^-) (11) 

the dynamical exponent z can be estimated by plot- 
ting (m(T)'')/(m(T)^)^ versus ri^^ at the known = 
2.2691 • • • value (Fig. The estimate of z is obtained 
when curves of different sizes are collapsed into a sin- 
gle curve. The best fit for data collapse is obtained 
by choosing z — 2.143 ± 0.063 for all the data shown 
in Fig. ^ For the curves of L = 32 and 48, the best 
fit occurs at z 2.11 ± 0.02. The best fit is ob- 
tained at z = 2.138 ± 0.039 for L = 48 and 64, and 



at z = 2.143 ± 0.063 for L = 64 and 80. We have 
given the error bars on z from the variance of the es- 
timate of z by a set of independent runs. Simulations 
were performed at relatively small system sizes compared 
to more extensive simulations [171 QM ■ Considering the 
effect of corrections to scaling, our value of z is, within 
error bars, approaching the values that were previously 
reported z = 2.16 - 2.17 17, 18, 19]. 

To summarize, we have presented a reweighting 
method for nonequilibrium Markov processes. We have 
shown the basis of this method starting from the for- 
mulation of the SIS method. As a demonstration, we 
have used the Ising model on a square lattice. With the 
nonequilibrium simulation at a single temperature, we 
can determine the critical temperature and critical expo- 
nents using finite-size scaling. The nonequilibrium simu- 
lation for large enough systems 0,0] is one way to extract 
dynamical and also static properties without worrying 
about finite-size effects. The systematic analysis of finite- 
size effects, the finite-size scaling is another 
way of studying nonequilibrium relaxation. The nonequi- 
librium reweighting is useful when it is combined with the 
latter approach. We should note that this method is nei- 
ther restricted to single spin flip updates nor the Ising 
model. This method is applicable, in principle, when- 
ever the ratio of probabifities P/3'{a-l_^_i\af)/Pp{a(^i\af) 
can be calculated numerically. It should be applicable 
with other Monte Carlo update schemes such as the clus- 
ter update schemes j23, as well as the A^-fold way 
method '25| . We should mention that Dickman "2^ pro- 
posed a similar reweighting method for the application 
to the contact process. However, our formalism is more 
extensive and general. 
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